Modelling phytoplankton-virus interactions: phytoplankton blooms and lytic virus transmission

A dynamic reaction–diffusion model of four variables is proposed to describe the spread of lytic viruses among phytoplankton in a poorly mixed aquatic environment. The basic ecological reproductive index for phytoplankton invasion and the basic reproduction number for virus transmission are derived to characterize the phytoplankton growth and virus transmission dynamics. The theoretical and numerical results from the model show that the spread of lytic viruses effectively controls phytoplankton blooms. This validates the observations and experimental results of Emiliana huxleyi-lytic virus interactions. The studies also indicate that the lytic virus transmission cannot occur in a low-light or oligotrophic aquatic environment.


Introduction
Phytoplankton are the world's most important aquatic producers.They play an essential role in biogeochemical cycles and strongly influence the abundance and biodiversity of aquatic communities.Light and nutrients are two essential resources for phytoplankton growth (Huisman et al. 2006;Klausmeier and Litchman 2001;Wang et al. 2007;Yoshiyama et al. 2009;Zhang et al. 2021).Phytoplankton absorb light energy to synthesize carbon dioxide and water into organic matter and release oxygen (Chen et al. 2015;Davies and Wang 2021).At the same time, they ingest various nutrients from the surrounding environment to preserve normal physiological metabolism (Vasconcelos et al. 2016;Zhang et al. 2018Zhang et al. , 2021)).These processes are important for achieving energy conversion and elemental cycles in nature and maintaining the carbon-oxygen balance of the atmosphere.
Viruses in aquatic ecosystems are generally small particles.The components of viruses are mainly nucleic acids and protein coats (Fuhrman 1999).This means that viruses can trigger the biosynthesis of the viral genome and protein only when they are parasitic in the host cells (Beretta and Kuang 1998;Edwards and Steward 2018;Fuhrman 1999;Fuhrman et al. 2011).Viruses are extremely abundant and widely distributed over oceans, lakes, and rivers (Suttle 2005).It is shown that viruses are the major pathogens and important causes of mortality for most aquatic organisms (Demory et al. 2021;Fuhrman 1999).As a result, viruses directly affect the structure and stability of aquatic communities.
It is recognized that viruses infect a significant proportion of phytoplankton and it is a major cause of the loss of phytoplankton (Demory et al. 2021;Kuhlisch et al. 2021;Suttle et al. 1990).According to the mechanism of virus transmission among phytoplankton, lytic viruses are considered to be one of the most common viruses (Beretta and Kuang 1998;Edwards and Steward 2018;Fuhrman 1999).The process of the lytic virus infection can be divided into the following steps.First, the virus contacts and adsorbs on phytoplankton cells by random diffusion, and then injects its nucleic acid into the cells.Second, the virus takes over the synthesis machinery of phytoplankton cells and produces viral genome and protein biosynthesis, which is needed for the viral offspring.Third, after the new virus is assembled, the lytic process ends with the lysis of the phytoplankton cell membrane, and then the virus particles are released into aquatic environments.This infection process indicates that lytic viruses destroy phytoplankton cells and reduce the biomass of phytoplankton.In view of the interrelationship between phytoplankton and lytic viruses, it is of great interest to explore the spread of lytic viruses among phytoplankton.
In this study, we propose a mathematical model to describe the spread of lytic viruses among phytoplankton.Here the aquatic environment under consideration is a poorly mixed water column.This suggests that both phytoplankton and virus distributions are spatially heterogeneous (Huisman et al. 2006;Klausmeier and Litchman 2001).Light comes from the water surface and nutrients come from the bottom of the water column (Ryabov et al. 2010;Yoshiyama et al. 2009;Zhang et al. 2021).Phytoplankton growth requires light and nutrients.Lytic viruses move randomly with turbulence and use phytoplankton cells as hosts to replicate and reproduce, eventually releasing a large number of lytic viruses when the cells rupture (Demory et al. 2021;Edwards and Steward 2018;Fuhrman 1999;Kuhlisch et al. 2021).One principal objective of the present paper is to model and elucidate the mechanism of lytic virus transmission among phytoplankton in a heterogeneous environment.
Several mathematical models have been introduced to investigate lytic virus transmission (Béchette et al. 2013;Beretta and Kuang 1998;Demory et al. 2021;Edwards and Steward 2018;Fuhrman et al. 2011).In these existing studies, the basic assumptions are that both phytoplankton and viruses are spatially uniformly distributed.
However, there is growing evidence that only shallow aquatic environments and epilimnion are well-mixed, while most aquatic ecosystems are poorly mixed.This means that phytoplankton generally exhibit strong spatial heterogeneity (Huisman et al. 2006;Klausmeier and Litchman 2001;Ryabov et al. 2010;Yoshiyama et al. 2009;Zhang et al. 2021).Our model consists of four dynamic reaction-diffusion equations.Its contribution is the ability to describe the spatially heterogeneous distribution of phytoplankton and viruses.This model also incorporates the effects of light and phytoplankton sinking relative to existing models of lytic viruses.Thus the reaction-diffusion model contains advection terms and a nonlocal structure.This increases the complexity of the model structure and the difficulty of model analysis.We will rigorously derive the basic ecological reproductive index for phytoplankton invasion and the basic reproduction numbers for lytic virus transmission by analyzing nonnegative steady states and some basic properties of solutions of the model.
Phytoplankton blooms are an important phenomenon in which phytoplankton biomass increases rapidly and significantly over a period of time (Chen et al. 2015;Kuhlisch et al. 2021).It has adverse effects on aquatic ecological environments.Phytoplankton blooms lead to poor water quality of aquatic resources, cause the death of aquatic organisms, and even threaten the health and safety of humans (Ho et al. 2019).It has become apparent that lytic viruses can influence phytoplankton biomass abundance from observations and some experiments (Demory et al. 2021;Edwards and Steward 2018;Fuhrman et al. 2011;Kuhlisch et al. 2021;Suttle et al. 1990).For example, Emiliana huxleyi is distributed worldwide and frequently forms large and dense blooms.These blooms are often terminated by the lytic virus transmission (Kuhlisch et al. 2021).Experiments have shown that about 50% of Emiliana huxleyi cells are infected by a large double-stranded DNA virus during blooms, and 25-100% Emiliana huxleyi deaths are related to the lytic virus infection (Fuhrman 1999;Kuhlisch et al. 2021).Another objective of this study is to examine these observations and experimental results theoretically and reveal the evolution trend in phytoplankton biomass and free lytic virus density with varying ecological factors based on the mathematical model described above.
The structure of the paper is organized as follows.In the next section, a mathematical model consisting of four reaction-diffusion equations is formulated to describe phytoplankton-virus interactions.By using the principal eigenvalue theory, bifurcation theory, and persistence theory, we analyze some basic properties of dynamic solutions and nonnegative steady states of the model in Sect.3. Two important basic indices are derived.In Sect.4, numerical bifurcation and time series diagrams are made to explore the role of lytic virus transmission for phytoplankton blooms, and the changes in phytoplankton biomass and free lytic virus density for varying environmental factors.An overview of the main conclusions and future research questions are presented in the last section.

Model formulation
The water depth coordinate is x and the time scale is t.Consider a poorly mixed water column with x = 0 at the water surface and x = x l at the bottom of the water column.The model consists of four reaction-diffusion equations describing the change of the concentrations of susceptible phytoplankton (S(x, t)), virus-infected phytoplankton (I (x, t)), free lytic virus particles (V (x, t)), and dissolved nutrients (N (x, t)).Their interaction relationship is shown in Fig. 1.The biological significance of variables and parameters in the model is summarized in Table 1.
According to the research work in Beretta and Kuang (1998); Edwards and Steward (2018); Fuhrman et al. (2011), we have the following assumptions: (A 1 ) Only susceptible phytoplankton can reproduce through photosynthesis and consumption of nutrients; (A 2 ) Virus-infected phytoplankton are removed by lysis before reproducing.The latency period is T from the infection to the lysis.The lytic virus reproduces inside phytoplankton cells during this period T ; (A 3 ) The lysis of virus-infected phytoplankton releases massive amounts of lytic virus particles.
Let S(x, t) and I (x, t) be the biomass density of susceptible phytoplankton and virus-infected phytoplankton, respectively.They have two different forms of movement in the water column.One is random movement in the vertical direction by turbulence with a diffusivity d p (Huisman et al. 2006;Klausmeier and Litchman 2001;Ryabov et al. 2010;Yoshiyama et al. 2009;Zhang et al. 2021).The other is directional movement including sinking or buoyant due to gravity or seeking more light with a velocity ω (Grover 2017;Klausmeier and Litchman 2001;Ryabov et al. 2010;Yoshiyama et al. 2009).The whole water column is a closed environment for phytoplankton.This means that S(x, t) and I (x, t) satisfy no-flux boundary conditions at endpoints x = 0 and x = x l .
Whether virus-infected phytoplankton can consume resources is still a disputed subject (Gourley and Kuang 2004).Here we introduce a parameter θ ∈ [0, 1], which represents the proportion of infected phytoplankton capable of resource consumption (Smith and Thieme 2012).By Assumptions (A 1 ) and (A 2 ), the growth of susceptible phytoplankton depends on light and nutrients.It is expressed as  Here the light intensity L(x, S + θ I ) following the Lambert-Beer law (Huisman et al. 2002) is given by since light is absorbed by water and phytoplankton above the point x.The reduction of susceptible phytoplankton biomass includes three parts: −μ s S (death and grazing Jäger et al. 2010;Vasconcelos et al. 2016;Wang et al. 2007;Zhang et al. 2021), −η(S + θ I )S (competition among phytoplankton for other growth resources such as inorganic carbon for photosynthesis Davies and Wang 2021;Hsu et al. 2017;Nie et al. 2016;Zhang et al. 2021), and −bβ SV (lytic virus infection into virus-infected phytoplankton (Beretta and Kuang 1998;Edwards and Steward 2018;Fuhrman et al. 2011)).The increase in virus-infected phytoplankton biomass comes from bβ SV and the loss is owing to the lysis with the lytic death rate −δV where δ = 1/T and T is the latency period from the infection to the lysis.Let V (x, t) denote the density of free lytic viruses in the water column.By Assumption (A 3 ), its increase is from the lysis release of infected phytoplankton with the lytic virus replication factor q, also known as "burst size" (Beretta and Kuang 1998;Edwards and Steward 2018;Fuhrman et al. 2011).The reduction in V is caused by death with a rate −μ v V and infectious consumption with a rate −β SV .The free lytic virus particles move randomly in the water column under the influence of turbulence with a diffusion rate d v .Neumann boundary conditions at x = 0 and x = x l mean that no free lytic virus particles enter or leave the water column.
The function N (x, t) describes dissolved nutrient concentration in the water column.The nutrient supply is through the nutrient exchange at the bottom of the water column (x = x l ) with a fixed nutrient input concentration N 0 and an exchange rate α (Klausmeier and Litchman 2001;Ryabov et al. 2010;Yoshiyama et al. 2009;Zhang et al. 2018).The nutrients are transported by turbulence in the water column with a diffusion rate d n .A no-flux boundary condition is imposed at x = 0 since there is no nutrient exchange at the water surface.The nutrient consumption in the water column by susceptible and virus-infected phytoplankton is described by a rate The above assumptions and analysis yield the following complete model for phytoplankton-virus interactions with light and nutrients 123 for x ∈ (0, x l ) and t > 0 with the boundary conditions and the initial conditions (2.3) Here ω ∈ R, θ ∈ [0, 1], and remaining parameters are assumed to be positive constants.Model (2.1) is a system of four reaction-diffusion-advection equations with a nonlocal structure.To explore the phytoplankton-virus interactions, we rigorously analyze the dynamic properties of model (2.1) including basic properties and behavior of solutions and existence and stability of nonnegative steady states.

Model analysis
The main purpose of this section is to explore the theoretical results of model (2.1)-(2.3).Some basic properties of solutions of model (2.1)-(2.3)are given in Sect.3.1.The study of nonnegative steady state solutions is in Sect.3.2.The dynamic numerical simulations are performed to explain and supplement our theoretical results in Sect.3.3.

Basic properties of solutions
Let X = C([0, x l ], R 4 ) denote the Banach space of all continuous functions defined on [0, x l ] with values in R 4 and the norm being the supremum norm.The feasible domain W for (2.1)-(2.3) is the positive cone in X : (3.1) Theorem 3.1 The system (2.1)-(2.3)possesses a unique classical solution in W for all t > 0 and it is dissipative.

Proof
The local existence and uniqueness of nonnegative classical solutions of system (2.1)-(2.3)follow from standard arguments (see Martin and Smith 1990).To obtain the global existence of the solutions, we only need to prove that the solutions of system (2.1)-(2.3)are dissipative.
From the N -equation in (2.1), we obtain From the S-equation, we have Adding the S-equation and the I -equation gives for x ∈ (0, x l ) and t > t 1 with the boundary condition Applying the parabolic comparison theorem, we get lim sup (3.4) For the > 0 above, there exists a t 2 > 0 such that for any t ≥ t 2 .From the V -equation in (2.1), we have Hence, Combining (3.2)-(3.5)shows that the solutions of system (2.1)-(2.3)are dissipative.This completes the proof.

Steady states
In order to investigate the spread of lytic viruses among phytoplankton, we explore nonnegative steady states of model (2.1)-(2.3).There are three types of steady state solutions of model (2.1)-( 2.3) as follows.
(ii) Disease-free steady state E 2 = (S 2 (x), 0, 0, N 2 (x)), where S 2 (x) and N 2 (x) satisfy on (0, x l ) with the boundary conditions (3.8) To characterize the dynamic behavior of model (2.1)-( 2.3), we define the basic ecological reproductive index for phytoplankton invasion.For h ∈ L ∞ ([0, x l ]), let λ 1 (d p , ω, x l , h(x)) denote the principal eigenvalue of (3.9)By Proposition 3.1 in Wang et al. (2019), λ 1 exists and it is unique.Moreover, The basic ecological reproductive index for phytoplankton invasion is defined as Note that μ * p is related to the stability of the unique extinction steady state E 1 (see (3.11)).The index R p measures the reproductive capacity of phytoplankton, and it describes the average number of new phytoplankton cells produced by per cubic meter of phytoplankton in a life cycle.
Proof It is clear that E 1 ≡ (0, 0, 0, N 0 ) exists uniquely.The local stability of E 1 is determined by the eigenvalue problem with the boundary conditions (3.12d) One can observe that λ is an eigenvalue of (3.11) if and only if λ is an eigenvalue of one of the following four operators and it has at least one eigenvalue greater than 0 if R p > 1.The above analysis shows that E 1 is locally asymptotically stable when R p < 1, and E 1 is unstable when R p > 1.
To obtain the global stability of E 1 when R p < 1, we only need to prove that it is globally attractive.For any > 0, from (3.2) and (3.5), there exists a t Let ξ be the first component of the positive eigenfunction of (3.11) corresponding to μ = μ * p satisfying S(x, t * 1 ) ≤ cξ(x) for some c > 0. Let Ŝ = Se −(ω/d p )x and ξ = ξ e −(ω/d p )x .It follows from the S-equation in model (2.1) that By the comparison theorem of parabolic systems, we have This shows that lim sup and is sufficiently small.Hence we can find a t since is sufficiently small.Similarly, we can also obtain Following Theorem 1.8 in Mischaikow et al. (1995) or Theorem 4.1 in Thieme (1992), the N -equation in (2.1) becomes for sufficiently large t.Thus and then E 1 is globally attractive.
Remark 3.4 1.The basic ecological reproductive index for phytoplankton invasion R p measures the viability of phytoplankton.R p < 1 means that phytoplankton go extinct and nutrients are evenly distributed in the water column.μ p = μ * p is a critical loss rate that determines whether phytoplankton can invade an aquatic ecosystem.2. From the structure of (3.10), the basic ecological reproductive index R p depends on important ecological parameters such as spatial factors d p , ω, the water column depth x l , nutrient concentration at the bottom of the water column N 0 and light intensity at the water surface L 0 .According to Theorems 3.2 and 3.4 in Hsu and Lou (2010), we conclude that d R p /dω < 0 and d R p /dx l < 0. By the monotonicity of λ 1 with respect to rg(N 0 ) f (L(x, 0)), d R p /d N 0 > 0 and d R p /d L 0 > 0. Thus, high N 0 and L 0 are conducive to phytoplankton invasion, while high ω and x l prevent phytoplankton invasion.The dependence of R p on d p is complicated.Both high and low d p may be detrimental to the survival of phytoplankton (see Huisman et al. 2002).
When R p > 1, the extinction state E 1 is unstable and phytoplankton can persist.To establish the existence of E 2 which has a positive phytoplankton mass, we first derive a priori estimate for nonnegative solutions of (3.6).
From the strong maximum principle and Hopf boundary lemma, we have and then N 2 (x l ) < N 0 .By the N -equation in (3.6), we obtain This implies that N 2 (x) > 0 on [0, x l ] from the maximum principle.Note that Combining its boundary conditions give N 2 (x) > 0 on (0, x l ).Thus, 0 From the S-equation in (3.6) and S 2 > 0, we get By applying the similar arguments of Lemma 2.2 in Pang et al. (2019), we can conclude that S 2 (x) ≤ ((rg We now prove the existence of disease-free steady state E 2 by showing E 2 bifurcating from the line of extinction state {(μ p , E 1 ) : 3) has at least one disease-free steady state E 2 .(ii) There exists an ε > 0 such that for each given s ∈ (0, ε) the bifurcating solution (μ p (s), S 2 (s, •), N 2 (s, •)) is locally asymptotically stable with respect to the susceptible phytoplankton-nutrient model (3. 16) The first part of the proof is divided into two steps.The first one is to obtain the existence of local bifurcation of E 2 .Applying the Crandall-Rabinowitz bifurcation theorem (see Theorem 1.7 in Crandall and Rabinowitz 1971), we show that there is a positive solution branch The second is to explore the global bifurcation structure of E 2 .That is to show that the disease-free steady state E 2 exists for all μ p ∈ (0, μ * p ) by using Theorem 3.3 and Remark 3.4 in Shi and Wang (2009).
In Theorem 3.6, we obtain the existence of E 2 for all 0 < μ p < μ * p , and the uniqueness and stability of E 2 is unknown.Numerical simulations suggest that the disease-free steady state E 2 is unique if R p > 1.
Next we show that when R p > 1 and some additional conditions are satisfied, solutions of model (2.1)-(2.3)are uniformly persistent and there exists at least one endemic steady state E 3 .To obtain the conclusions, we define the basic reproduction number for lytic virus transmission.Letting Ī = I e −(ω/d p )x and linearizing (2.1) at a disease-free steady state E 2 = (S 2 , 0, 0, N 2 ), we obtain For h 1 , h 2 ∈ C([0, x l ], R + ), we consider the following linear parabolic system R 2 ) denote the solution semigroup generated by the following system We assume that the distribution of initial infected phytoplankton and free lytic viruses is ρ(x) = (ρ 1 (x), ρ 2 (x)).As time evolves, the distribution at time t is h 2 (t)ρ(x).Therefore, it can be deduced that the distribution of total new infected phytoplankton is Here K h 1 ,h 2 is called the next generation operator, and its spectral radius is r (K h 1 ,h 2 ).It follows that the basic reproduction number associated with (h 1 , h 2 ) is given as Especially, if h 1 = h 2 = S 2 , then the basic reproduction number associated with E 2 for virus transmission is denoted as Since the uniqueness of E 2 is not known, R 0 defined in (3.22) depends on E 2 , and the definition given in (3.21) allows a more general basic reproduction number which will be used in the following.
Applying the similar arguments of Theorem 3.1 (i) in Wang and Zhao (2012), one can obtain the following conclusion.
From Theorem 3.3, we have lim t→∞ (S(•, t), N (•, t)) = (0, N 0 ) when R p < 1.The numerical simulations in the next section suggest that the attractor 0 when R p > 1 only contains the steady state (S 2 , N 2 ).Let and Here W is defined in (3.1).For any (S, N ) ∈ U, we introduce a projection H by (3.25) The uniform persistence shown in Lemma 3.8 implies that 0 Let ϒ 0 := {(S, 0, 0, N ) ∈ W : (S, N ) ∈ 0 }.It will prove that E 1 and ϒ 0 are uniform weak repellers with respect to W * .Lemma 3.9 If R p > 1, then E 1 is a uniform weak repeller for (2.1)-( 2.3) with respect to W * , that is, there is a υ 1 > 0 satisfying lim sup t→∞ dist( (t)σ 0 , E 1 ) ≥ υ 1 for any Note that both f and g are continuous.Thus, for the above > 0, there is a υ 1 > 0 such that Assume that the conclusion is not true.Then we can find a σ 0 ∈ W * satisfying lim sup (3.27)This suggests that for some sufficiently large T 1 .It follows from the strong maximum principle and the Hopf boundary lemma that S(•, T 1 , σ 0 ) > 0 since σ 0 ∈ W * .Hence, there exists a x and ξ = ξ e −(ω/d p )x .From (3.26) and the S-equation in (2.1), we have It is easy to see that From the comparison theorem of parabolic system, we obtain It contradicts (3.27).Therefore, E 1 is a uniform weak repeller with respect to W * .
Since 0 is compact, we can find an St ∈ 0 satisfying Hence, From the I -and V -equations in (2.1), we let Ī = I e −(ω/d p )x and have By the strong maximum principle and the Hopf boundary lemma, Ī (x, T 2 , σ 0 ) > 0 and Applying the comparison theorem of parabolic system again, we have This implies that Ī (x, t, σ 0 ), V (x, t, σ 0 ) are unbounded since λ 0 (S * − υ 2 , S * + υ 2 ) > 0. This contradicts with (3.28).The proof is completed.Now we are ready to prove the existence of E 3 and uniform persistence of model (2.1)-(2.3)when R p > 1 and R 0 (S * , S * ) > 1.
Theorem 3.11 Let S * and S * be defined as in (3.25) 3) possesses at least one endemic steady state E 3 .
From numerical simulations, we speculate that the global attractor 0 = {(S 2 , N 2 )} in Lemma 3.8 and the basic reproduction number R 0 (S * , S * ) = R 0 (S 2 , S 2 ) = R 0 can be uniquely determined.Theorem 3.11 shows that lytic viruses will be transmitted in the phytoplankton population if R p > 1 and R 0 > 1. Hence R 0 = 1 is a threshold for lytic viruses from persistence to extinction. Figure 2 reveals the evolution of R 0 with respect to model parameters.From numerical simulations in Sect.3.3, the detailed lytic virus transmission is more complicated.The endemic state may be a positive steady state E 3 , or a positive spatially inhomogeneous periodic solution.In Fig. 3, the extinct steady state E 1 is the attractor with R p = 0.49.Phytoplankton are extinct and thus the lytic virus does not spread in the water column.Nutrients tend to the input concentration N 0 and are spatially uniformly distributed.Figure 4 shows that the (possibly unique) disease-free steady state E 2 is the attractor when R p = 4.91 and R 0 = 0.13.Phytoplankton invade aquatic ecosystems and exhibit vertical aggregation.Lytic virus still cannot be transmitted among phytoplankton because lytic virus mortality is high.
The model (2.1)-(2.3) is uniformly persistent and an endemic steady state E 3 appears to be the asymptotic state when R p = 4.91 and R 0 = 1.59 (see Fig. 5).As μ v decreases further (R p = 4.91, R 0 = 3.24), a positive spatially inhomogeneous periodic solution becomes the asymptotic state, which arises from bifurcate from E 3 via a Hopf bifurcation (see Fig. 6).Both Figs. 5 and 6 show the lytic virus is prevalent in the phytoplankton, either in a form of spatially heterogenous steady state (Fig. 5) or a spatially heterogenous temporal-oscillatory fashion (Fig. 6).It is noticeable that susceptible phytoplankton, virus-infected phytoplankton and lytic viruses aggregate near the water surface (for only short time in the oscillatory case).Comparing Fig. 4 (a 2 )-(b 2 ) and Fig. 5 (a 3 )-(b 3 ), one can also observe that lytic viruses reduce phytoplankton biomass.

Phytoplankton blooms and lytic virus transmission
Phytoplankton blooms are an important manifestation of the pollution of the aquatic environments, and even lead to the collapse of entire aquatic ecosystems.It has been shown that the wide spreading of lytic viruses transmission among phytoplankton can control phytoplankton blooms from observations and experiments (Fuhrman 1999;Kuhlisch et al. 2021).In the following discussion, we will verify these observations and experimental results with the proposed model (2.1)-(2.3).It is noted that the role of ecological factors in the interaction of phytoplankton blooms and lytic virus transmission is not very clear.It is also necessary and significant to explore the effects of ecological factors in this process.
We focus on the environmental parameters related to viral transmission and ecological factors in model (2.1)-(2.3).Those parameters include viral infection-related rates β, q, spatial ecological factors d p , d v , d n , ω, and resource-related ecological factors L 0 and N 0 .In the figures below, we show the evolution of asymptotic states (steady state solutions E 1 , E 2 , E 3 and the spatially inhomogeneous periodic solution) for different parameter values.When the solutions of (2.1)-(2.3)converge to one of E 1 , E 2 , E 3 or the spatially inhomogeneous periodic solution, numerical bifurcation diagrams show the evolution trend of densities of susceptible phytoplankton ((1/x l ) x l 0 S(x)dx), virus-infected phytoplankton ((1/x l ) x l 0 I (x)dx), and free lytic viruses ((1/x l ) x l 0 V (x)dx).For the time-periodic solutions, the minimum and maximum values are shown.Time series diagrams reveal the evolution of densities of susceptible phytoplankton (((1/x l ) x l 0 S(x, t)dx), infected phytoplankton ((1/x l ) x l 0 I (x, t)dx) and lytic viruses ((1/x l ) x l 0 V (x, t)dx) over time.The parameter values of the numerical analysis used here are derived from Table 1.The simulations are implemented in MATLAB via the finite difference method.
We first examine the effect of lytic virus transmission parameters β, q.Changes in these parameters are closely related to ecological factors such as temperature, salinity (Demory et al. 2021).The transmission coefficient β is an important indicator for  2.3) with the transmission coefficient β.For β = 0, one can observe that lytic viruses do not spread among phytoplankton and that susceptible phytoplankton biomass is at a high level.When β increases, susceptible phytoplankton biomass gradually declines, while lytic virus loads ascend sharply and a significant proportion of phytoplankton are infected by viruses.In Fig. 7a, b, spatially inhomogeneous periodic solutions bifurcate from positive steady states through a Hopf bifurcation at β = 0.0032.This shows that increasing the transmission coefficient will lead to persistent phytoplankton blooms, but a large transmission coefficient will cause a large amplitude pulse bloom which occurs periodically about every thirty days.In the latter case, an increase of susceptible phytoplankton precedes the peak of virus transmission, then the phytoplankton population is at a low level for a long time until the next wave.
The lytic virus burst size varies with host genotype or environmental conditions.It is an important factor for describing the transmission of lytic viruses among phytoplankton (Edwards and Steward 2018).Figure 8 shows an evolution of asymptotic states in model (2.1)-(2.3)for varying burst size q similar to the one for transmission coefficient.The dynamics of the model progress from a disease-free steady state, to an endemic steady state, and then to a spatially inhomogeneous periodic solution.This indicates that the small burst size does not cause phytoplankton infection and that only a certain degree of burst size can lead to the spread of lytic viruses among phytoplankton and a decrease in total phytoplankton biomass.
Emiliana huxleyi are eukaryotic microscopic phytoplankton that are widely distributed in oceans or brackish waters.They are one of the important producers in marine ecosystems, especially playing an important role in the carbon cycle of oceans.Emiliana huxleyi have frequent, large-scale blooms each year.It is known that Emiliana huxleyi can be massively infected by a lytic virus with a large double-stranded DNA structure.The lytic virus is called Emiliana huxleyi virus and causes massive mortality of Emiliana huxleyi.In Fuhrman (1999); Kuhlisch et al. (2021), the authors stated that Emiliana huxleyi blooms are often terminated by Emiliana huxleyi virus from experiments and observations.The periodic pattern shown in Fig. 7b is consistent with the observations and experimental results of Emiliana huxleyi-lytic virus interactions in Fuhrman (1999); Kuhlisch et al. (2021).At the beginning of a bloom cycle, both phytoplankton and lytic virus densities are at a very low value.After that, the susceptible phytoplankton biomass rises rapidly and reaches a maximum, while the virus density remains almost constant.Phytoplankton blooms occur at this time.As the time progresses further, there is a sudden and dramatic increase in lytic virus density.This indicates a large-scale spread of lytic viruses among phytoplankton.During this process, the total phytoplankton biomass declines sharply.When the virus reaches its maximum, phytoplankton biomass tends to almost zero, thus the bloom ends.In the last stage of this bloom cycle, the viral density decreases, and is again at a low value together with phytoplankton. Figure 8b also exhibits a similar periodic bloom cycle pattern.The above findings effectively validate the large-scale transmission of lytic viruses to terminate phytoplankton blooms.
The evolution of phytoplankton and lytic virus densities with different spatial ecological factors (turbulent diffusivity and directional movement velocity) can be seen in Figs. 9 and 10.Phytoplankton, viruses, and nutrients all move randomly in the water column with turbulence, so it is assumed that they have the same diffusion coefficient If there are no viruses in the aquatic environment, the phytoplank-Fig.10 Effect of the sinking or buoyant velocity ω on the density of susceptible phytoplankton, virusinfected phytoplankton and lytic viruses Fig. 11 Effects of the water surface light intensity L 0 and sediment input nutrient concentration N 0 on the density of susceptible phytoplankton, virus-infected phytoplankton and lytic viruses ton biomass rises rapidly and reaches a high value with increasing turbulence intensity.It is because adequate nutrient transport in the water column facilitates better phytoplankton growth (see Fig. 9a).When the lytic virus spreads among phytoplankton, the biomass of phytoplankton, including susceptible and infected phytoplankton, remains almost unchanged for D ∈ (0.02, 5) except the periodic oscillatory patterns occurring for the intermediate diffusion rate (see Fig. 9b).In this process, there are two stability switches from steady states to periodic oscillations, and then back to steady states.This illustrates that the presence of lytic viruses not only leads to complex dynamics, but also effectively reduces the phytoplankton biomass.Comparing Fig. 10a, b, similar conclusions are obtained for the directional movement velocity ω.This once again confirms that lytic virus transmission can effectively control phytoplankton blooms.
Light and nutrients are two important ecological factors, and their abundance directly affects the invasion of phytoplankton and the spread of lytic viruses.We choose the water surface light intensity L 0 and sediment nutrient input concentration N 0 as typical resource-related ecological factors to explore the survival and extinction of phytoplankton and lytic viruses.Figure 11 shows that low L 0 or N 0 causes phytoplankton extinction since light and nutrients are two essential resources for phy-Fig.12 Parameter regions of L 0 versus N 0 for the survival and extinction of phytoplankton and lytic viruses.Phytoplankton are extinct in 1 (R p < 1), phytoplankton survive while lytic viruses become extinct in 2 (R p > 1, R 0 < 1), lytic viruses spread among phytoplankton in a steady state form in 3 (R p > 1, R 0 > 1) or as a spatially inhomogeneous periodic solution form in 4 (R p > 1, R 0 > 1) toplankton growth.With the gradual increase of L 0 or N 0 , phytoplankton first invade aquatic ecosystems, after which lytic viruses begin to spread among phytoplankton.When L 0 or N 0 is at a high value, the system undergoes a Hopf bifurcation and exhibits periodic oscillations.This is a classical paradox of enrichment.
The combined effects of L 0 and N 0 are presented in Fig. 12.Using the basic ecological reproductive index R p and basic reproduction number R 0 to distinguish dynamic behavior, we take L 0 and N 0 as the coordinate parameters to divide the first quadrant of the L 0 − N 0 plane into four regions i for i = 1, 2, 3, 4. Phytoplankton go extinct in 1 (R p < 1); Phytoplankton successfully invade aquatic habitats, while lytic viruses cannot survive in 2 (R p > 1 and R 0 < 1).If R p > 1 and R 0 > 1, lytic viruses spread among phytoplankton as a steady state form in 3 or as a spatially inhomogeneous periodic solution form in 4 .These findings suggest that it is very difficult for lytic viruses to spread in a low-light or oligotrophic aquatic ecosystem.

Conclusion and discussion
Viruses are pervasive components of aquatic ecosystems, and have important influences on aquatic biodiversity and biogeochemical cycles (Suttle 2005;Suttle et al. 1990).As one of the common viruses, lytic viruses use phytoplankton cells as their primary hosts and are ubiquitous in aquatic ecosystems (Demory et al. 2021;Edwards and Steward 2018;Fuhrman et al. 2011;Kuhlisch et al. 2021).They replicate and reproduce inside phytoplankton cells, and eventually release virions by lysing the cells.The spread of lytic viruses causes a dramatic decrease in phytoplankton biomass and even terminates phytoplankton blooms.Few mathematical models have been formulated to explore the interactions between phytoplankton and lytic viruses and to expound the mechanisms of lytic virus transmission.
The dynamic model (2.1)-(2.3) is proposed to describe the spread of lytic virus among phytoplankton in a water column.Compared with the existing models in Béchette et al. (2013); Beretta and Kuang (1998); Demory et al. (2021); Edwards and Steward (2018); Fuhrman et al. (2011), there are two novel points in model (2.1)-(2.3).One is to consider a poorly mixed aquatic environment and take account of random movements of phytoplankton and viruses.The other is to add phytoplankton sinking movement and include light as a factor contributing to the growth of phytoplankton.Our results show that these newly added factors have a significant effect on the dynamics of phytoplankton growth and the spread of lytic viruses.
Two quantitative ecological indices: the basic ecological reproductive index R p for phytoplankton invasion and the basic reproduction number R 0 for virus transmission, are rigorously derived from the model.According to Theorems 3.3, 3.6, 3.11 and corresponding remarks, phytoplankton go extinct if R p < 1, phytoplankton successfully invade and lytic viruses are extinct if R p > 1 and R 0 < 1, lytic viruses spread among phytoplankton as a steady state form or a spatially inhomogeneous periodic solution form if R p > 1 and R 0 > 1.By using theoretical and numerical analysis of the model (2.1)-( 2.3), we consider the interaction between phytoplankton blooms and lytic virus transmission and the role of ecological factors.From the numerical bifurcation diagrams (Figs. 7,8,9,10,11), one can observe that the spread of lytic viruses controls phytoplankton biomass.Time series diagrams (Figs. 7, 8b) reveal large-scale outbreaks of lytic virus infection effectively terminate phytoplankton blooms.The findings validate the observations and experimental results of Emiliana huxleyi-lytic virus interactions in Fuhrman (1999); Kuhlisch et al. (2021).The studies also indicate that it is very difficult for lytic viruses to spread in a low-light or oligotrophic aquatic environment (see Figs. 11,12).
In this paper, we attempt to model the spread of lytic virus among phytoplankton, and the mechanism of lytic virus transmission and its relationship with phytoplankton blooms are explored.There are still many mathematical and biological problems worthy of further study.Mathematically, more dynamic properties of model (2.1)-(2.3)need to be rigorously investigated.For example, the uniqueness and stability of disease-free steady state E 2 and endemic steady state E 3 , and the existence of spatially inhomogeneous periodic solutions.Biologically, the phytoplankton intracellular nutrient-to-carbon ratio in model (2.1)-(2.3) is assumed to be constant, but in reality, it is significantly variable (Loladze et al. 2000;Wang et al. 2007).It is desirable to include this factor in phytoplankton-virus interactions.Bacteria are a very important aquatic microorganism.They decompose organic carbon produced by phytoplankton and are infected by aquatic viruses (Fuhrman 1999;Suttle 2005;Wang et al. 2007;Yan et al. 2022).It is very natural to add bacteria into the model (2.1)-(2.3)and further explore the carbon cycle in aquatic ecosystems.In addition, the effects of zooplankton, fish (Chen et al. 2017) and toxins (Shan and Huang 2019;Nie et al. 2017) could also be considered.

Fig. 1
Fig. 1 Phytoplankton-virus interactions in a poorly mixed water column S t = d p S xx − ωS x diffusion and advection + rg(N ) f (L(x, S + θ I ))S susceptible phytoplankton growth − μ p S − η(S + θ I )S virus-infected phytoplankton , V t = d v V xx diffusion + qδ I lysis release of virus-infected phytoplankton − μ v V free lytic virus death − β SV infection , N t = d n N xx diffusion − c p rg(N ) f (L(x, S + θ I ))(S + θ I ) 18) and the first equation in (3.19), we have x l 0 κ 1 (x)e (ω/d p )x ξ(x)dx = 0, and ξ ∈ W 1 can be uniquely solved under this condition.It follows from the Fredholm alternative theorem that ζ ∈W 2 can then be uniquely solved by the second equation in (3.19).Hence, range in (3.20) and (3.23) respectively, we can define the principal eigenvalue λ 0 (S * , S * ) and the basic reproduction number R 0 (S * , S * ) associated with (S * , S * ).It follows from Lemma 3.7 that λ 0 (S * , S * ) and R 0 (S * , S * ) − 1 have the same sign.

Fig. 7 a
Fig. 7 a Effect of the transmission coefficient β on the density of susceptible phytoplankton, virus-infected phytoplankton and lytic viruses.b Time series of density of phytoplankton and viruses

Fig. 8 aFig. 9
Fig. 8 a Effect of the burst size q on the density of susceptible phytoplankton, virus-infected phytoplankton and lytic viruses.b Time series of density of phytoplankton and viruses

Table 1
Variables and parameters with realistic values and biological significance of model